#ifndef CPARTICLE_H
#define CPARTICLE_H

#include <vector>
#include <iostream>
#include <stdio.h>
#include <string>
#include <math.h>
#include <time.h>
#include <list>
using namespace std;

#include "../Point/CPoint.h"
#include "../Heap/Heap.h"
#include "../CLineScan/CLineScan.h"

#define PI 3.14159265359

template <typename _TO>
class CPattern {
	vector<_TO>		_features;
	
public:
	string _label;

	CPattern (int _nf = 0){
		for (int i = 0; i < _nf; i++)
			_features.push_back(0.0);
	}

	void push_back(const _TO& _f){
		_features.push_back(_f);
	}

	void operator += (const CPattern& _p){
		for (int i = 0; i < (int)_p._features.size(); i++)
			_features[i] += _p._features[i];
	}

	double operator () (void){
		double _average = 0.0;
		for (int i = 0; i < (int)_features.size(); i++)
			_average +=  _features[i];
		_average /= _features.size();
		return _average;
	}

	bool operator <= (const CPattern& _p){
		bool _less = true;
		for (int i = 0; i < (int)_p._features.size(); i++)
			if (_features[i] > _p._features[i])
				_less = false;
		return _less;
/*
		double _norm = NULL, __norm = NULL;
		for (int i = 0; i < (int)_p._features.size(); i++){
			_norm += _features[i] * _features[i];
			__norm += _p._features[i] * _p._features[i];
		}
		return sqrt(_norm) <= sqrt(__norm); 
*/	}

	bool operator < (const CPattern& _p){
		double _norm = NULL, __norm = NULL;
		for (int i = 0; i < (int)_p._features.size(); i++){
			_norm += _features[i] * _features[i];
			__norm += _p._features[i] * _p._features[i];
		}
		return sqrt(_norm) < sqrt(__norm); 
	}


	bool operator == (const CPattern& _p){
		bool _equal = true;
		for (int i = 0; i < (int)_p._features.size(); i++)
			if (_p._features[i] != _features[i]) _equal =  false;
		return _equal;
	}

	CPattern operator + (double _f){
		CPattern _p;
		for (int i = 0; i < (int)_p._features.size(); i++)
			_p.push_back(_features[i] + _f);
		return _p;
	}

	CPattern& operator + (const CPattern& _p){
		for (int i = 0; i < (int)_p._features.size(); i++)
			_features[i] += _p._features[i];
		return *this;
	}
				

	CPattern& operator * (int _f){
		for (int i = 0; i < (int)_features.size(); i++)
			_features[i] *= _f;
		return *this;
	}

	CPattern& operator / (int _f){
		for (int i = 0; i < (int)_features.size(); i++)
			_features[i] /= _f;
		return *this;
	}

	long double operator - (const CPattern& _p){
		long double dist = 0.0;
		for (int i = 0; i < (int)_p._features.size(); i++)
			dist += (_features[i] - _p._features[i])*(_features[i] - _p._features[i]);
		return sqrt(dist);
	}

	double norma (void){
		double _norm = NULL;
		for (int i = 0; i < (int)_features.size(); i++){
			_norm += _features[i] * _features[i];
		}
		return sqrt(_norm); 
	}

	int size (void){
		return (int)_features.size();
	}

	void operator /= (int _d){
		for (int i = 0; i < (int)_features.size(); i++)
			_features[i] /= _d;
	}
	
	template <typename __TO>
	friend  ostream& operator << (ostream& out, const CPattern <__TO>& p){
		out << "[ ";
		for (int i = 0; i  < (int)p._features.size(); i++)
			out << p._features[i] << " ";
		out << "]";
		return out;
	}
};

class CParticleTraits {
public:
	typedef double					value_pattern_t;
	typedef int						value_point_t;
	typedef CPoint<value_point_t>	point_t;
	typedef CPattern<value_pattern_t>	pattern_t;
	typedef double					gravity_t;
	typedef double					velocity_t;
	typedef double					angle_t;
	typedef string					label_t;
};


template <typename _CT> class CParticle{
public:
	typedef typename _CT::pattern_t			pattern_t;
	typedef typename _CT::point_t			point_t;
	typedef typename _CT::gravity_t			gravity_t;
	typedef typename _CT::velocity_t		velocity_t;
	typedef typename _CT::angle_t			angle_t;
	typedef typename _CT::value_pattern_t	value_pattern_t;
	typedef typename _CT::value_point_t		value_point_t;
	typedef typename _CT::label_t			label_t;
	typedef CLineScan<CParticle,point_t>	line_scan_t;

	typedef list<pattern_t>				container_t;
//	typedef Min<CParticle>					min_t;
//	typedef Heap<CParticle,min_t>			container_t;
	
public:
	pattern_t	_pattern;
	point_t		_point_in_plane;
	container_t	_container; // particle container y masa
	velocity_t	_velocity;
	gravity_t	_gravitational_focus; // radio gravitacional
	angle_t		_angle;
	label_t		_label;
	bool		_is_planetesimal;
	bool		_locked;

	CParticle ***	_nebula;  // Asociated Nebula
	int				_limit;
	double			_velocity_factor;
	double			_gravity_factor;
	pattern_t		_mean;
	double			_mean_error;
	double			_covariance;

public:
	CParticle (const velocity_t& _v = NULL, const pattern_t& _p = pattern_t())
	:	_velocity(_v), _angle(NULL), _pattern(_p), _limit(NULL),
		_gravitational_focus(NULL), _label("NULL"), _is_planetesimal(false), _locked(false){
	}


	CParticle (const CParticle& p)
	: _limit(NULL), _angle(NULL){
		_velocity = p._velocity;
		_gravitational_focus = p._gravitational_focus;
		_is_planetesimal = p._is_planetesimal;
		_locked = p._locked;
		_container = p._container;
	}

	~CParticle (void) {
	}

	double operator() (void){
		return _pattern();
	}


	void update_values (void){
		// Setting  mean
		_mean = pattern_t(_mean.size());
		container_t::iterator it = _container.begin();
		for (; it != _container.end(); it++)
			_mean += *it;
		_mean /= _container.size();
		
//		_covariance = 0.0;
//		it = _container.begin();
//		for (; it != _container.end(); it++)
//			_covariance += _mean - 	*it;
//		_covariance /= _container.size();
	}

	void operator() (point_t& _p){ // Scan  Line Functor
		if (_point_in_plane == _p) return;
		if (!_nebula[_point_in_plane[X]][_point_in_plane[Y]]) return;
		if (_nebula[_p[X]][_p[Y]]){
			if (_nebula[_p[X]][_p[Y]]->_is_planetesimal == false){
//				if (_nebula[_p[X]][_p[Y]]->_pattern <= _mean){
				if (_nebula[_p[X]][_p[Y]]->_pattern <= _mean + _mean_error){
//				cout << _nebula[_p[X]][_p[Y]]->_pattern - _mean << " " <<  _covariance << endl;
//				if (_nebula[_p[X]][_p[Y]]->_pattern - _mean <= _covariance){
					(*_nebula[_point_in_plane[X]][_point_in_plane[Y]]) << _nebula[_p[X]][_p[Y]]->_pattern;
					if (_velocity > 0) 	_velocity -= _velocity_factor;
					_gravitational_focus += _gravity_factor;
					_mean = ((_mean * (_container.size() - 1)) 
							+ _nebula[_p[X]][_p[Y]]->_pattern) 
							/ _container.size();
					
//					update_values();
					delete _nebula[_p[X]][_p[Y]];
					_nebula[_p[X]][_p[Y]] = NULL;

				}				
			}
		}
//		update_values();

		int range	= (int)floor(_gravitational_focus + .5);
		int init_x	= (_p[X] - range > 0) ? _p[X] - range : 0;
		int fin_x	= (_p[X] + range >= _limit) ? _limit - 1 : _p[X] + range;
		int init_y	= (_p[Y] - range > 0) ? _p[Y] - range : 0;
		int fin_y	= (_p[Y] + range >= _limit) ? _limit - 1 : _p[Y] + range;

		for (int i = init_x; i < fin_x + 1; i++){
			for (int k = init_y; k < fin_y + 1; k++){
				if (_nebula[i][k]){				
					if (_nebula[i][k]->_is_planetesimal == true) continue;
					if (_point_in_plane == point_t(i,k)) continue;
//					if (_nebula[i][k]->_pattern <= _mean){
					if (_nebula[i][k]->_pattern <= _mean + _mean_error){
//					if (_nebula[i][k]->_pattern - _mean <= _covariance){
						(*_nebula[_point_in_plane[X]][_point_in_plane[Y]]) << _nebula[i][k]->_pattern;
						if (_velocity > 0) 	_velocity -= _velocity_factor;
						_gravitational_focus += _gravity_factor;
						_mean = ((_mean * (_container.size() - 1)) + _nebula[i][k]->_pattern) / _container.size();
						delete _nebula[i][k];
						_nebula[i][k] = NULL;
//						update_values();
					}
				}
			}
		}
//		update_values();
	}

	void acretion (CParticle ***& _n, int _l, double _gf, double _vf, double _merr){
		_gravity_factor = _gf;
		_velocity_factor = _vf;
		_mean_error = _merr;
		_limit = _l;
		_nebula = _n;
		if (_velocity > 0){
			_angle = PI * (rand() % 360) / 180;		
			int nx = (int)floor((_velocity * cos(_angle)) + .5);
			int ny = (int)floor((_velocity * sin(_angle)) + .5);
			nx = (nx < 0) ? 0 : nx;
			ny = (ny < 0) ? 0 : ny;
			nx = (nx >= _limit) ? _limit - 1 : nx;
			ny = (ny >= _limit) ? _limit - 1 : ny;
			
			while (_nebula[nx][ny]){
				if (_nebula[nx][ny]->_is_planetesimal == false) break;
				_angle = PI * (rand() % 360) / 180;
				nx = (int)floor((_velocity * cos(_angle)) + .5);
				ny = (int)floor((_velocity * sin(_angle)) + .5);
				nx = (nx < 0) ? 0 : nx;
				ny = (ny < 0) ? 0 : ny;
				nx = (nx >= _limit) ? _limit - 1 : nx;
				ny = (ny >= _limit) ? _limit - 1 : ny;
			}
			line_scan_t (_point_in_plane,point_t(nx,ny),this);
			swap(_nebula[_point_in_plane[X]][_point_in_plane[Y]],_nebula[nx][ny]);				
			_point_in_plane[X] = nx;
			_point_in_plane[Y] = ny;
		}
	}


	void lock_pattern (void){
		_locked = true;
//		_container.push_back(this->_pattern);
		_mean = (this->_pattern);
//		_mean = pattern_t(this->_pattern.size());
		_covariance = _mean.norma();
	}

	CParticle& operator << (const CParticle& _p){
		_container.push_back(_p->_pattern);
		return *this;
	}

	CParticle& operator << (const pattern_t& _p){
		_container.push_back(_p);
		return *this;
	}


	CParticle& operator << (const value_pattern_t& _v){
		if (!_locked) _pattern.push_back(_v);
		return *this;
	}

	CParticle& operator << (const string& _l){
		_label = _l;
		this->_pattern._label = _l;
		return *this;
	}

	value_point_t& operator[](const int& pos){
		return _point_in_plane[pos];
	}

	template <typename __CT> 
	friend ostream& operator << (ostream& out, const CParticle<__CT>& _p){
		out << "Parcticle: ";
		out <<"P: " << _p._pattern;
		out << "\tV: " << _p._velocity;
		out << "\tF: " << _p._gravitational_focus;
		out << "\tM: " << _p._container.size();
		out << "\tI: " << _p._is_planetesimal;
		out << "\tP: " << _p._point_in_plane;
		return out;
	}

};

#endif